Sleep

0. Библиотеки и предварительная подготовка

library(readxl)
library(dplyr)
library(scatterPlotMatrix)
library(kableExtra)
library(psych)
library(summarytools)
library(ggplot2)
library(GGally)
library(nortest)
library(moments)
library(ppcor)
library(corrplot)
library(lm.beta)
library(ellipse)
library(car)
library(olsrr)
print_df <- function(df)
{
  df |>
    kable(format = "html") |>
    kable_styling() |>
    kableExtra::scroll_box(width = "100%", height = "100%")
}

plot_qq_graph <- function(data, column_name) {
  data<-data
  expected_quantiles <- qnorm(ppoints(length(data)))
  qqnorm(data,  pch = 19, col = "deeppink")
  qqline(data, distribution = qnorm, lwd = 2, col = "limegreen")
  legend("topleft", legend = column_name, col = "deeppink", pch = 19)
}

plot_pp_plot <- function(data, column_name) {
  n <- length(data)
  expect_prob <- pnorm(data, mean(data), sqrt(var(data) * (n - 1) / n))
  expect_prob <- sort(expect_prob)
  plot(x = expect_prob, y = ppoints(n), col = "red")
  lines(x = c(0, 1), y = c(0, 1), col = "blue")
  legend("topleft", legend = column_name, col = "deeppink", pch = 19)
}

1. Загрузка и предобработка данных

df <- read_excel("stat2021_sm/Sleep/SLEEP_shortname.xls", )
df <- df |> mutate(SLEEP = ifelse(SLEEP < 0, NA, SLEEP), PARADOX = ifelse(is.infinite(PARADOX), NA, PARADOX), PRED_IND = as.factor(PRED_IND), EXP_IND = as.factor(EXP_IND), DANG_IND = as.factor(DANG_IND)) 
df |> slice(1:10) |> print_df()

NAME	BODY_WEI	BRAIN_WE	SLOWWAVE	PARADOX	SLEEP	LIFESPAN	GESTTIME	PRED_IND	EXP_IND	DANG_IND
African elephant	6654.000	5712.0	NA	NA	3.3	38.6	645	3	5	3
African giant rat	1.000	6.6	6.3	2.0	8.3	4.5	42	3	1	3
Arctic Fox	3.385	44.5	NA	NA	12.5	14.0	60	1	1	1
Arc. ground squirrel	0.920	5.7	NA	NA	16.5	NA	25	5	2	3
Asian elephant	2547.000	4603.0	2.1	1.8	3.9	69.0	624	3	5	4
Baboon	10.550	179.5	9.1	0.7	9.8	27.0	180	4	4	4
Big brown bat	0.023	0.3	15.8	3.9	19.7	19.0	35	1	1	1
Brazilian tapir	160.000	169.0	5.2	1.0	6.2	30.4	392	4	5	4
Cat	3.300	25.6	10.9	3.6	14.5	28.0	63	1	2	1
Chimpanzee	52.160	440.0	8.3	1.4	9.7	50.0	230	1	1	1

2. Описание признаков

1. NAME — название животного — качественный признак
2. BODY_WEI — вес тела в килограммах — количественны непрерывный признак
3. BRAIN_WE — вес мозга в граммах — количественный непрерывный признак
4. SLOWWAVE — время медленного сна, часов в день — количественный непрерывный признак
5. PARADOX —время быстрого сна, часов в день — количественный непрерывный признак
6. SLEEP — общее время сна, часов в день — количественный непрерывный признак
7. LIFESPAN — продолжительность жизни в годах — количественный непрерывный признак
8. GESTTIME — время беременности в днях — количественный непрерывный признак
9. PRED_IND — индекс хищничества (1 — наименее вероятный объект для охоты, 5 — наиболее вероятный объект для охоты) — порядковый признак
10. EXP_IND — индекс воздействия во время сна (1 — наименее подверженный воздействию (например, животное спит в хорошо защищенном логове), 5 — наиболее подверженный воздействию) — порядковый признак
11. DANG_IND — общий индекс опасности (на основе двух вышеуказанных индексов и другой информации) (1 — наименьшая опасность (от других животных), 5 — наибольшая опасность (со стороны других животных)) — порядковый признак

3. Матричный график разброса

categories <-list(NULL, NULL, NULL, NULL, NULL, NULL, NULL, 1:5, 1:5, 1:5)
df |> dplyr::select(-NAME) |> scatterPlotMatrix(regressionType = 1,
                        corrPlotType = "Text",
                        slidersPosition = list(
                          dimCount = 7,
                          xStartingDimIndex = 1,
                          yStartingDimIndex = 1
                        ),
                        categorical = categories,
                        plotProperties = list(noCatColor = "Indigo"),
                        controlWidgets = TRUE,
                        height = 1050, width = 1000)

Видим достаточно логичные кандидаты на логарифмирование: вес тела, вес мозга, продолжительность жизни и время беременности — это мультиплекативные признаки. Также время быстрого сна имеет перекос влево, так что лучше его тоже отлогарифмировать, хотя логически, так как время сна содержит в себе время быстрого сна, этот признак не является мультиплекативным, а также не имеет нелинейную зависимость. Однако, ассиметрия времи быстрого сна имеет достаточно большое значение:

desc0 <- df |>dplyr::select(-NAME, -PRED_IND, -EXP_IND, -DANG_IND) |> describe() |>dplyr::select(-trimmed, -mad, -se)
print_df(desc0)

	vars	n	mean	sd	median	min	max	range	skew	kurtosis
BODY_WEI	1	62	198.789984	899.158011	3.3425	0.005	6654.0	6653.995	6.2494291	40.5981819
BRAIN_WE	2	62	283.134193	930.278942	17.2500	0.140	5712.0	5711.860	4.8288287	23.2367752
SLOWWAVE	3	48	8.672917	3.666452	8.3500	2.100	17.9	15.800	0.2798669	-0.4442636
PARADOX	4	50	1.972000	1.442651	1.8000	0.000	6.6	6.600	1.3667228	1.7793264
SLEEP	5	58	10.532759	4.606760	10.4500	2.600	19.9	17.300	0.1909556	-0.6497094
LIFESPAN	6	58	19.877586	18.206255	15.1000	2.000	100.0	98.000	1.9107968	5.0043782
GESTTIME	7	58	142.353448	146.805039	79.0000	12.000	645.0	633.000	1.5974241	2.3223297

4. Логарифмирование хвостатых

df.log <- df
df.log <- df.log |> mutate(BODY_WEI = log(BODY_WEI), BRAIN_WE = log(BRAIN_WE), PARADOX = log(PARADOX), LIFESPAN = log(LIFESPAN), GESTTIME = log(GESTTIME))
df.log <- df.log |> mutate(PARADOX = ifelse(is.infinite(PARADOX), NA, PARADOX))
df.log |> dplyr::select(-NAME) |> filter(!is.na(SLOWWAVE) & !is.na(PARADOX) & !is.na(SLEEP) & !is.na(LIFESPAN) & !is.na(GESTTIME)) |> scatterPlotMatrix(regressionType = 1,
                        corrPlotType = "Text",
                        slidersPosition = list(
                          dimCount = 7,
                          xStartingDimIndex = 1,
                          yStartingDimIndex = 1
                        ),
                        categorical = categories,
                        plotProperties = list(noCatColor = "Indigo"),
                        controlWidgets = TRUE,
                        height = 1050, width = 1000)

Выведем также немного другой матричный график, в котором удаление NA будет попарным.

df.log |> dplyr::select(-NAME) |> ggpairs(columns = 1:7)

Видим, что распределения стали более симметричными, причём почти все признаки имеют отрицательный эксцесс, то есть они более “плоские”, чем нормальное:

desc1 <- df.log |>dplyr::select(-NAME, -PRED_IND, -EXP_IND, -DANG_IND) |> describe() |>dplyr::select(-trimmed, -mad, -se)
print_df(desc1)

	vars	n	mean	sd	median	min	max	range	skew	kurtosis
BODY_WEI	1	62	1.3375390	3.1231277	1.2066382	-5.2983174	8.802974	14.101291	0.1453599	-0.5232692
BRAIN_WE	2	62	3.1401979	2.4465120	2.8477071	-1.9661129	8.650324	10.616437	0.0434421	-0.6184785
SLOWWAVE	3	48	8.6729167	3.6664517	8.3500000	2.1000000	17.900000	15.800000	0.2798669	-0.4442636
PARADOX	4	49	0.4656547	0.7080779	0.5877867	-1.2039728	1.887070	3.091042	-0.1364307	-0.5840339
SLEEP	5	58	10.5327586	4.6067604	10.4500000	2.6000000	19.900000	17.300000	0.1909556	-0.6497094
LIFESPAN	6	58	2.5930586	0.9435145	2.7120343	0.6931472	4.605170	3.912023	-0.1629152	-0.9151912
GESTTIME	7	58	4.4441910	1.0602517	4.3596587	2.4849066	6.469250	3.984344	0.0473325	-1.1407247

По графику неравномерностей не видно. Если раскрасить по категориальным признакам, то тоже не получается выделить раздельные облака.

5. Выбросы

Видим пару выбросов, для которых вес мозга мал, но продолжительность жизни достаточно велика:

df.log |> filter(LIFESPAN > 2 & BRAIN_WE < 0) |> print_df()

NAME	BODY_WEI	BRAIN_WE	SLOWWAVE	PARADOX	SLEEP	LIFESPAN	GESTTIME	PRED_IND	EXP_IND	DANG_IND
Big brown bat	-3.772261	-1.203973	15.8	1.3609766	19.7	2.944439	3.555348	1	1	1
Little brown bat	-4.605170	-1.386294	17.9	0.6931472	19.9	3.178054	3.912023	1	1	1

Это две летучие мыши.

Построим теперь графики разброса без выбросов.

df.log.out <- df.log |>  filter(!(LIFESPAN > 2 & BRAIN_WE < 0))
df.log.out |> dplyr::select(-NAME) |> filter(!is.na(SLOWWAVE) & !is.na(PARADOX) & !is.na(SLEEP) & !is.na(LIFESPAN) & !is.na(GESTTIME)) |> scatterPlotMatrix(regressionType = 1,
                        corrPlotType = "Text",
                        slidersPosition = list(
                          dimCount = 7,
                          xStartingDimIndex = 1,
                          yStartingDimIndex = 1
                        ),
                        categorical = categories,
                        plotProperties = list(noCatColor = "Indigo"),
                        controlWidgets = TRUE,
                        height = 1050, width = 1000)

И сравним корреляции.

df.log |> dplyr::select(-NAME) |> ggpairs(columns = 1:7)

df.log.out |> dplyr::select(-NAME) |> ggpairs(columns = 1:7)

В основном они выросли, по крайней мере, там, где мы исключали выбросы.

6. Зависимости от индексов опасности

Необходимо описать разницу между животными по индексу опасности места, где они спят. Посмотрим на количество животных в каждой группе:

df.log.out |>dplyr::select(EXP_IND) |> summary()

 EXP_IND
 1:25   
 2:13   
 3: 4   
 4: 5   
 5:13   

Выборки для индексов \(2\) и \(5\) можем считать сбалансированными.

Визуально оценим распределения с помощью ящиков с усами:

df.log.out |> ggplot(aes(x = EXP_IND, y = BODY_WEI)) +
  geom_boxplot()

df.log.out |> ggplot(aes(x = EXP_IND, y = BRAIN_WE)) +
  geom_boxplot()

df.log.out |> ggplot(aes(x = EXP_IND, y = PARADOX)) +
  geom_boxplot()

df.log.out |> ggplot(aes(x = EXP_IND, y = SLOWWAVE)) +
  geom_boxplot()

df.log.out |> ggplot(aes(x = EXP_IND, y = SLEEP)) +
  geom_boxplot()

df.log.out |> ggplot(aes(x = EXP_IND, y = LIFESPAN)) +
  geom_boxplot()

df.log.out |> ggplot(aes(x = EXP_IND, y = GESTTIME)) +
  geom_boxplot()

В основном видна монотонная зависимость, однако, предстоит проверить это с помощью критериев.

7. Нормальность, вид распределения

Из-за малости выборок в третьей и четвёртой группах исключим их из исследования.

get_df_on_ind <- function(i) {
  return(df.log.out |> filter(EXP_IND == i))
}

Посмотрим сначала на гистограммы распределений по группам для каждого признака.

for (i in c(1, 2, 5)) {
  print(paste("Группа ", i))
  par(mfrow = c(2, 4))
  hist(get_df_on_ind(i)$BODY_WEI, main = "")
  hist(get_df_on_ind(i)$BRAIN_WE, main = "")
  hist(get_df_on_ind(i)$SLOWWAVE, main = "")
  hist(get_df_on_ind(i)$PARADOX, main = "")
  hist(get_df_on_ind(i)$SLEEP, main = "")
  hist(get_df_on_ind(i)$LIFESPAN, main = "")
  hist(get_df_on_ind(i)$GESTTIME, main = "")
}

[1] "Группа  1"
[1] "Группа  2"

[1] "Группа  5"

Посмотрим на нормальность с помощью Normal Probability Plot:

for (i in c(1, 2, 5)) {
  print(paste("Группа ", i))
  par(mfrow = c(2, 4))
  plot_qq_graph(get_df_on_ind(i)$BODY_WEI, "BODY_WEI")
  plot_qq_graph(get_df_on_ind(i)$BRAIN_WE, "BRAIN_WE")
  plot_qq_graph(get_df_on_ind(i)$SLOWWAVE, "SLOWWAVE")
  plot_qq_graph(get_df_on_ind(i)$PARADOX, "PARADOX")
  plot_qq_graph(get_df_on_ind(i)$SLEEP, "SLEEP")
  plot_qq_graph(get_df_on_ind(i)$LIFESPAN, "LIFESPAN")
  plot_qq_graph(get_df_on_ind(i)$GESTTIME, "GESTTIME")
}

[1] "Группа  1"
[1] "Группа  2"

[1] "Группа  5"

И с помощью PP-plot:

for (i in c(1, 2, 5)) {
  print(paste("Группа ", i))
  par(mfrow = c(2, 4))
  plot_pp_plot((get_df_on_ind(i) |> filter(!is.na(BODY_WEI)))$BODY_WEI, "BODY_WEI")
  plot_pp_plot((get_df_on_ind(i) |> filter(!is.na(BRAIN_WE)))$BRAIN_WE, "BRAIN_WE")
  plot_pp_plot((get_df_on_ind(i) |> filter(!is.na(SLOWWAVE)))$SLOWWAVE, "SLOWWAVE")
  plot_pp_plot((get_df_on_ind(i) |> filter(!is.na(PARADOX)))$PARADOX, "PARADOX")
  plot_pp_plot((get_df_on_ind(i) |> filter(!is.na(SLEEP)))$SLEEP, "SLEEP")
  plot_pp_plot((get_df_on_ind(i) |> filter(!is.na(LIFESPAN)))$LIFESPAN, "LIFESPAN")
  plot_pp_plot((get_df_on_ind(i) |> filter(!is.na(GESTTIME)))$GESTTIME, "GESTTIME")
}

[1] "Группа  1"
[1] "Группа  2"

[1] "Группа  5"

А теперь оценим по тесту Шапиро-Уилка:

df.shapiro.test.p.value <- data.frame(ind = c(1, 2, 5), BODY_WEI_P = rep(0, 3), BRAIN_WE_P = rep(0, 3), SLOWWAVE_P = rep(0, 3), PARADOX_P = rep(0, 3), SLEEP_P = rep(0, 3), LIFESPAN_P = rep(0, 3), GESTTIME_P = rep(0, 3))

for (i in c(1, 2, 5)) {
  df.shapiro.test.p.value[df.shapiro.test.p.value$ind == i, ] <- c(i, shapiro.test(get_df_on_ind(i)$BODY_WEI)$p.value, 
  shapiro.test(get_df_on_ind(i)$BRAIN_WE)$p.value,
  shapiro.test(get_df_on_ind(i)$SLOWWAVE)$p.value,
  shapiro.test(get_df_on_ind(i)$PARADOX)$p.value,
  shapiro.test(get_df_on_ind(i)$SLEEP)$p.value,
  shapiro.test(get_df_on_ind(i)$LIFESPAN)$p.value,
  shapiro.test(get_df_on_ind(i)$GESTTIME)$p.value)
}

df.shapiro.test.p.value |> print_df()

ind	BODY_WEI_P	BRAIN_WE_P	SLOWWAVE_P	PARADOX_P	SLEEP_P	LIFESPAN_P	GESTTIME_P
1	0.1924485	0.4855703	0.4685556	0.0946346	0.3176217	0.0733372	0.2755650
2	0.0204661	0.0126712	0.4403360	0.8802017	0.4562102	0.8049845	0.5314537
5	0.9922347	0.7461442	0.0844404	0.9214902	0.0052604	0.4522227	0.0232529

И наконец по Лиллиефорсу:

df.lillie.test.p.value <- data.frame(ind = c(1, 2, 5), BODY_WEI_P = rep(0, 3), BRAIN_WE_P = rep(0, 3), SLOWWAVE_P = rep(0, 3), PARADOX_P = rep(0, 3), SLEEP_P = rep(0, 3), LIFESPAN_P = rep(0, 3), GESTTIME_P = rep(0, 3))

for (i in c(1, 2, 5)) {
  df.lillie.test.p.value[df.lillie.test.p.value$ind == i, ] <- c(i, lillie.test(get_df_on_ind(i)$BODY_WEI)$p.value, 
  lillie.test(get_df_on_ind(i)$BRAIN_WE)$p.value,
  lillie.test(get_df_on_ind(i)$SLOWWAVE)$p.value,
  lillie.test(get_df_on_ind(i)$PARADOX)$p.value,
  lillie.test(get_df_on_ind(i)$SLEEP)$p.value,
  lillie.test(get_df_on_ind(i)$LIFESPAN)$p.value,
  lillie.test(get_df_on_ind(i)$GESTTIME)$p.value)
}

df.lillie.test.p.value |> print_df()

ind	BODY_WEI_P	BRAIN_WE_P	SLOWWAVE_P	PARADOX_P	SLEEP_P	LIFESPAN_P	GESTTIME_P
1	0.3977150	0.3724944	0.2399545	0.2094751	0.2837660	0.3215462	0.4647484
2	0.1209793	0.0760938	0.8789409	0.9154588	0.8195576	0.7890638	0.9238088
5	0.7872499	0.4314403	0.0328885	0.5556462	0.0004679	0.4954892	0.1272532

Этот критерий является менее мощным, однако даёт схожие результаты.

По итогу можно считать, что распределения всех признаков во всех группах похожи на нормальные, кроме веса тела и веса мозга во второй группе и времени медленного сна, времени сна и времени беременности в пятой группе.

8. t-тесты

Сначала посмотрим на дисперсию, чтобы применить t-test для распределений с равными дисперсиями. Тест Фишера проверяет равенство дисперсий для нормально распределённых данных, которые мы определелили в прошлом пункте:

df.fisher.var <- data.frame(IND1 = rep(0, 3), IND2 = rep(0, 3), BODY_WEI_P = rep(0, 3), BRAIN_WE_P = rep(0, 3), SLOWWAVE_P = rep(0, 3), PARADOX_P = rep(0, 3), SLEEP_P = rep(0, 3), LIFESPAN_P = rep(0, 3), GESTTIME_P = rep(0, 3))

k <- 1

for (i in c(1, 2, 5)){
  for (j in c(1, 2, 5)){
    if (i < j){
      df.fisher.var[k, ] <- c(i, j,
                                  var.test(get_df_on_ind(i)$BODY_WEI, get_df_on_ind(j)$BODY_WEI)$p.value,
                                  var.test(get_df_on_ind(i)$BRAIN_WE, get_df_on_ind(j)$BRAIN_WE)$p.value,
                                  var.test(get_df_on_ind(i)$SLOWWAVE, get_df_on_ind(j)$SLOWWAVE)$p.value,
                                  var.test(get_df_on_ind(i)$PARADOX, get_df_on_ind(j)$PARADOX)$p.value,
                                  var.test(get_df_on_ind(i)$SLEEP, get_df_on_ind(j)$SLEEP)$p.value,
                                  var.test(get_df_on_ind(i)$LIFESPAN, get_df_on_ind(j)$LIFESPAN)$p.value,
                                  var.test(get_df_on_ind(i)$GESTTIME, get_df_on_ind(j)$GESTTIME)$p.value)
      k <- k + 1
    }
  }
}

df.fisher.var |> print_df()

IND1	IND2	BODY_WEI_P	BRAIN_WE_P	SLOWWAVE_P	PARADOX_P	SLEEP_P	LIFESPAN_P	GESTTIME_P
1	2	0.3119679	0.1905823	0.3041400	0.4869923	0.6645465	0.9543801	0.7543563
1	5	0.6380604	0.3940484	0.4168502	0.7837348	0.0334406	0.0080450	0.6722005
2	5	0.6237712	0.6698967	0.1563733	0.4297079	0.0221523	0.0125095	0.5226789

Гипотеза о равенстве не отвергается для веса тела и веса мозга между первой и пятой группы, для времени медленного сна, времени сна и времени беременности между первой и второй группы, времени быстрого сна между всеми группами и продолжительности жизни между первой и второй группой.

Проведём t-тесты для первой, второй и пятой группы по всем признакам:

df.t.test.p.value <- data.frame(IND1 = rep(0, 3), IND2 = rep(0, 3), BODY_WEI_P = rep(0, 3), BRAIN_WE_P = rep(0, 3), SLOWWAVE_P = rep(0, 3), PARADOX_P = rep(0, 3), SLEEP_P = rep(0, 3), LIFESPAN_P = rep(0, 3), GESTTIME_P = rep(0, 3))

k <- 1
flag <- FALSE
flagBB <- FALSE
flagSSG <-FALSE
flagL <- FALSE

for (i in c(1, 2, 5)){
  for (j in c(1, 2, 5)){
    if (i < j){
        flag <- i == 2 & j == 5
        flagBB <- i == 1 & j == 5
        flagSSG <- i == 1 & j == 2
        flagL <- i == 1 & j == 2
      
      df.t.test.p.value[k, ] <- c(i, j,
                                  t.test(get_df_on_ind(i)$BODY_WEI, get_df_on_ind(j)$BODY_WEI, var.equal = flag | flagBB)$p.value,
                                  t.test(get_df_on_ind(i)$BRAIN_WE, get_df_on_ind(j)$BRAIN_WE, var.equal = flag | flagBB)$p.value,
                                  t.test(get_df_on_ind(i)$SLOWWAVE, get_df_on_ind(j)$SLOWWAVE, var.equal = flag | flagSSG)$p.value,
                                  t.test(get_df_on_ind(i)$PARADOX, get_df_on_ind(j)$PARADOX, var.equal = TRUE)$p.value,
                                  t.test(get_df_on_ind(i)$SLEEP, get_df_on_ind(j)$SLEEP, var.equal = flag | flagSSG)$p.value,
                                  t.test(get_df_on_ind(i)$LIFESPAN, get_df_on_ind(j)$LIFESPAN, var.equal = flag | flagL)$p.value,
                                  t.test(get_df_on_ind(i)$GESTTIME, get_df_on_ind(j)$GESTTIME, var.equal = flag | flagSSG)$p.value)
      k <- k + 1
    }
  }
}

df.t.test.p.value |> print_df()

IND1	IND2	BODY_WEI_P	BRAIN_WE_P	SLOWWAVE_P	PARADOX_P	SLEEP_P	LIFESPAN_P	GESTTIME_P
1	2	0.3331031	0.8379613	0.4313905	0.0894578	0.3325090	0.4306579	0.2064574
1	5	0.0000011	0.0000094	0.0000195	0.0000105	0.0000000	0.0000080	0.0000807
2	5	0.0000004	0.0000050	0.0045885	0.0168625	0.0000441	0.0026227	0.0033338

Гипотеза о равенстве распределений по всем признакам между первой и второй не отвергается, а между первой и пятой и второй и пятой отвергается. Однако для признаков, которые мы посчитали ненормальными (вес тела и вес мозга во второй группе и время медленного сна, время сна и время беременности в пятой группе) p-value нельзя считать верным, так как размеры выборок малы и эти признаки не имеют нормального распределения.

9. Тест Манна-Уитни

Для проверки равенства распределений для ненормально распределённых признаков можем применить критерий Манна-Уитни, однако прежде посмотрим на симметричност по группам:

df.skew <- data.frame(IND = c(1, 2, 5), BODY_WEI_P = rep(0, 3), BRAIN_WE_P = rep(0, 3), SLOWWAVE_P = rep(0, 3), PARADOX_P = rep(0, 3), SLEEP_P = rep(0, 3), LIFESPAN_P = rep(0, 3), GESTTIME_P = rep(0, 3))

for (i in c(1, 2, 5)){
  df.skew[df.skew$IND == i, ] <- c(i, skewness(get_df_on_ind(i)$BODY_WEI, na.rm = TRUE),
                                      skewness(get_df_on_ind(i)$BRAIN_WE, na.rm = TRUE),
                                      skewness(get_df_on_ind(i)$SLOWWAVE, na.rm = TRUE),
                                      skewness(get_df_on_ind(i)$PARADOX, na.rm = TRUE),
                                      skewness(get_df_on_ind(i)$SLEEP, na.rm = TRUE),
                                      skewness(get_df_on_ind(i)$LIFESPAN, na.rm = TRUE),
                                      skewness(get_df_on_ind(i)$GESTTIME, na.rm = TRUE))
}

df.skew |> print_df()

IND	BODY_WEI_P	BRAIN_WE_P	SLOWWAVE_P	PARADOX_P	SLEEP_P	LIFESPAN_P	GESTTIME_P
1	0.1672298	0.5645030	-0.0002586	-0.5598674	0.3874961	1.0385890	0.5038765
2	-1.4951275	-1.6495435	0.5921447	-0.1315432	0.1049490	0.0122452	0.1394734
5	-0.1198309	0.0708243	1.0734624	-0.1861544	1.5538884	0.5578426	-1.0319607

Видим, что асимметрия веса тела и веса мозга у второй группы (ненормально распределённые признаки) сильно (примерно на порядок) отличается от первой и пятой группы. Аналогично: асимметрия времени сна и времени беременности пятой группы сильно отличается от первой и пятой. Для асимметрии времени быстрого сна у пятой группы и первой также нет оснований говорить о приблизительном равенстве, однако между пятой и второй группой разница мала, меньше \(0.5\), что может говорить об удовлетворении условий для критерия, однако надо ещё посмотреть на исправленную выборочную дисперсию:

df.var <- data.frame(IND = c(1, 2, 5), BODY_WEI_P = rep(0, 3), BRAIN_WE_P = rep(0, 3), SLOWWAVE_P = rep(0, 3), PARADOX_P = rep(0, 3), SLEEP_P = rep(0, 3), LIFESPAN_P = rep(0, 3), GESTTIME_P = rep(0, 3))

for (i in c(1, 2, 5)){
  df.var[df.var$IND == i, ] <- c(i, var(get_df_on_ind(i)$BODY_WEI, na.rm = TRUE),
                                    var(get_df_on_ind(i)$BRAIN_WE, na.rm = TRUE),
                                    var(get_df_on_ind(i)$SLOWWAVE, na.rm = TRUE),
                                    var(get_df_on_ind(i)$PARADOX, na.rm = TRUE),
                                    var(get_df_on_ind(i)$SLEEP, na.rm = TRUE),
                                    var(get_df_on_ind(i)$LIFESPAN, na.rm = TRUE),
                                    var(get_df_on_ind(i)$GESTTIME, na.rm = TRUE))
}

df.var |> print_df()

IND	BODY_WEI_P	BRAIN_WE_P	SLOWWAVE_P	PARADOX_P	SLEEP_P	LIFESPAN_P	GESTTIME_P
1	6.034299	4.526511	7.279905	0.3159571	12.914130	0.9033046	0.8095168
2	3.443109	2.189395	12.449889	0.4525376	15.658077	0.9042128	0.6633179
5	4.599834	2.815790	3.733333	0.2533922	3.156556	0.1911274	0.9828249

Видим, что дисперсии времени медленного сна для пятой и второй группы даже приблизительно не равны, поэтому ни один из признаков, которые мы хотели проверить не подходят для применения критерия Манна-Уитни. То есть для них мы ничего не можем сказать о равенстве распределений.

10. Коэффициент Пирсона

pairwise.cor <- function(data, name_method = "pearson"){
  m <- cor(data)
  m.p.value <- cor(data)
  
  for (i in colnames(data)){
    for (j in colnames(data)){
      data.out <- data|> filter(!is.na(data[[i]]) & !is.na(data[[j]]))
      m[i, j] <- cor(x = data.out[[i]], y = data.out[[j]], method = name_method)
      m.p.value[i, j] <- cor.test(x = as.vector(data.out[[i]]), y = as.vector(data.out[[j]]), method = name_method)$p.value
    }
  }
  
  corrplot(m, method = "number")
  as.data.frame(m.p.value)
}

Перед вычислением коэффициента корреляции Пирсона посмотрим на нормальность признаков для всех индивидов по тесту Шапиро-Уилка:

df.shapiro.test.all.p.value <- data.frame(BODY_WEI_P = 0, BRAIN_WE_P = 0, SLOWWAVE_P = 0, PARADOX_P = 0, SLEEP_P = 0, LIFESPAN_P = 0, GESTTIME_P = 0)

for (i in c(1)) {
  df.shapiro.test.all.p.value[1, ] <- c(shapiro.test(df.log.out$BODY_WEI)$p.value, 
  shapiro.test(df.log.out$BRAIN_WE)$p.value,
  shapiro.test(df.log.out$SLOWWAVE)$p.value,
  shapiro.test(df.log.out$PARADOX)$p.value,
  shapiro.test(df.log.out$SLEEP)$p.value,
  shapiro.test(df.log.out$LIFESPAN)$p.value,
  shapiro.test(df.log.out$GESTTIME)$p.value)
}

df.shapiro.test.all.p.value |> print_df()

BODY_WEI_P	BRAIN_WE_P	SLOWWAVE_P	PARADOX_P	SLEEP_P	LIFESPAN_P	GESTTIME_P
0.6745281	0.7912255	0.7502073	0.8786641	0.2600607	0.334608	0.1265908

Как видим по p-value можем считать, что распределения похожи на нормальные, поэтому критерий проверки на значимость коэффициента корреляции будет проверять ещё и независимость признаков, а также коэффициент Спирмена будет примерно равен коэффициенту Пирсона.

Теперь посчитаем коэффициенты корреляции Пирсона для всех пар признаков, причём удаление NA будем производить попарно, чтобы для каждой пары получить максимальное количество индивидов:

df.log.out |> 
  mutate(BRAIN_BODY = log(exp(BRAIN_WE) / exp(BODY_WEI))) |>
  dplyr::select(-NAME, -EXP_IND, -PRED_IND, -DANG_IND) |> 
  pairwise.cor() |> 
  print_df()

	BODY_WEI	BRAIN_WE	SLOWWAVE	PARADOX	SLEEP	LIFESPAN	GESTTIME	BRAIN_BODY
BODY_WEI	0.0000000	0.0000000	0.0003011	0.0182618	0.0002898	0.0000000	0.0000000	0.0000000
BRAIN_WE	0.0000000	0.0000000	0.0004412	0.0043862	0.0000914	0.0000000	0.0000000	0.0000397
SLOWWAVE	0.0003011	0.0004412	0.0000000	0.0000332	0.0000000	0.0008653	0.0000535	0.0414927
PARADOX	0.0182618	0.0043862	0.0000332	0.0000000	0.0000000	0.0060667	0.0000022	0.7322638
SLEEP	0.0002898	0.0000914	0.0000000	0.0000000	0.0000000	0.0000957	0.0000000	0.1188045
LIFESPAN	0.0000000	0.0000000	0.0008653	0.0060667	0.0000957	0.0000000	0.0000000	0.0241044
GESTTIME	0.0000000	0.0000000	0.0000535	0.0000022	0.0000000	0.0000000	0.0000000	0.0144867
BRAIN_BODY	0.0000000	0.0000397	0.0414927	0.7322638	0.1188045	0.0241044	0.0144867	0.0000000

Видим, что для всех пар признаков критерий даёт p-value меньше \(\alpha = 0.05\), то есть отвергаем гипотезу о незначимости зависимсоти для этих признаков.

Зависимости:

1. Вес мозга зависит от веса тела;

2. Общая продолжительность сна зависит от медленного и быстрого сна, а также они зависят друг от друга, но скорее косвенно, так как это соотношение скорее индивидуально для каждого животного;

3. Время беременности зависит от продолжительности жизни, так как беременность не должна занимать большей части жизни из-за уязвимости такого животного;

4. Чем больше животное, тем дольше беременность и продолжительность жизни, так как крупное животное требует как большого формирования во время беременности, так как и большего времени взросления;

5. Также от размера мозга зависит продолжительность жизни, однако эта зависимость не является прямой, так как если мы рассмотрим долю веса мозга к весу тела, то корреляция значительно упадёт;

6. Зависимости между сном и весом тела, сном и продолжительностью жизни или сном и времени беременности достаточно сложно объяснить, для этого нужно углублятся в сложные биологичесике процессы.

11. Коэффициент Спирмена

Проверим, что ранговый коэффициент корреляции Спирмена примерно равен Пирсону:

df.log.out |> 
 dplyr::select(-NAME, -EXP_IND, -PRED_IND, -DANG_IND) |> 
  pairwise.cor(name_method = "spearman") |> 
  print_df()

	BODY_WEI	BRAIN_WE	SLOWWAVE	PARADOX	SLEEP	LIFESPAN	GESTTIME
BODY_WEI	0.0000000	0.0000000	0.0019593	0.0063690	0.0004095	0.0000000	0.00e+00
BRAIN_WE	0.0000000	0.0000000	0.0008980	0.0004787	0.0000531	0.0000000	0.00e+00
SLOWWAVE	0.0019593	0.0008980	0.0000000	0.0000730	0.0000000	0.0010380	1.59e-05
PARADOX	0.0063690	0.0004787	0.0000730	0.0000000	0.0000000	0.0026538	2.60e-06
SLEEP	0.0004095	0.0000531	0.0000000	0.0000000	0.0000000	0.0000883	1.00e-07
LIFESPAN	0.0000000	0.0000000	0.0010380	0.0026538	0.0000883	0.0000000	0.00e+00
GESTTIME	0.0000000	0.0000000	0.0000159	0.0000026	0.0000001	0.0000000	0.00e+00

Результаты схожи и это не удивительно, ведь мы считаем, что данные похожи на нормальные, а для них отличие коэффициента Пирсона от Спирмена незначительно.

12. Частные корреляционные отношения

Однако возникает подозрение, что многие зависимости имеют другую причину, нежели зависимость только друг от друга. Как было видно по ящикам с усами для индекса опасности места для сна есть монотонная зависимость почти всех признаков от этого индекса. Поэтому возникает идея проверки коэффициента частной корреляции с вычетом влияния этого индекса.

df.log.out |> 
  group_by(EXP_IND) |> 
  mutate(BODY_WEI = BODY_WEI - mean(BODY_WEI, na.rm = TRUE), BRAIN_WE = BRAIN_WE - mean(BRAIN_WE, na.rm = TRUE), SLOWWAVE = SLOWWAVE - mean(SLOWWAVE, na.rm = TRUE), PARADOX = PARADOX - mean(PARADOX, na.rm = TRUE), SLEEP = SLEEP - mean(SLEEP, na.rm = TRUE), LIFESPAN = LIFESPAN - mean(LIFESPAN, na.rm = TRUE), GESTTIME = GESTTIME - mean(GESTTIME, na.rm = TRUE)) |>
  group_by(.drop = TRUE) |> dplyr::select(-NAME, -EXP_IND, -PRED_IND, -DANG_IND) |> 
  pairwise.cor() |> 
  print_df()

	BODY_WEI	BRAIN_WE	SLOWWAVE	PARADOX	SLEEP	LIFESPAN	GESTTIME
BODY_WEI	0.0000000	0.0000000	0.1044414	0.3927193	0.5207584	0.0000000	0.0000003
BRAIN_WE	0.0000000	0.0000000	0.0715185	0.9971036	0.1841235	0.0000000	0.0000000
SLOWWAVE	0.1044414	0.0715185	0.0000000	0.0055510	0.0000000	0.0723565	0.0053589
PARADOX	0.3927193	0.9971036	0.0055510	0.0000000	0.0000016	0.7878256	0.0094343
SLEEP	0.5207584	0.1841235	0.0000000	0.0000016	0.0000000	0.0607429	0.0002778
LIFESPAN	0.0000000	0.0000000	0.0723565	0.7878256	0.0607429	0.0000000	0.0000633
GESTTIME	0.0000003	0.0000000	0.0053589	0.0094343	0.0002778	0.0000633	0.0000000

Здесь мы видим, что достаточно много зависимостей перестали быть значимыми (\(\alpha = 0.05\)), так как мы исключили влияние индекса опасности места для сна. Ожидаемо зависимости всех видов сна от веса тела, веса мозга и продолжительности сна значительно ослабли, однако зависимости сна от времени беременности уменьшились мало.

13. Множественная регрессия

Посмотрим на данные:

df.log.out |> dplyr::select(-NAME) |> ggpairs(columns = 1:7)

Всё, что нужно, было прологарифмировано раньше, поэтому приступим к регрессии. Хотим прогнозировать время жизни по остальным данным.

Строим регрессию по всем признакам:

df.log.out.na <- df.log.out |> dplyr::filter(!is.na(SLOWWAVE) & !is.na(PARADOX) & !is.na(SLEEP) & !is.na(LIFESPAN) & !is.na(GESTTIME))
model <- lm(LIFESPAN ~ ., df.log.out.na |> dplyr::select(-NAME))
model.beta <- lm.beta(model)
summary(model.beta)

Call:
lm(formula = LIFESPAN ~ ., data = dplyr::select(df.log.out.na, 
    -NAME))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.77116 -0.16275  0.02196  0.12238  0.51269 

Coefficients:
             Estimate Standardized Std. Error t value Pr(>|t|)    
(Intercept)  1.348343           NA   0.903753   1.492 0.151321    
BODY_WEI    -0.022443    -0.065240   0.096158  -0.233 0.817825    
BRAIN_WE     0.461812     1.096994   0.104491   4.420 0.000264 ***
SLOWWAVE     0.095508     0.327255   0.176317   0.542 0.594021    
PARADOX      0.121882     0.080915   0.338311   0.360 0.722426    
SLEEP       -0.050191    -0.215930   0.181925  -0.276 0.785466    
GESTTIME    -0.099790    -0.105443   0.159882  -0.624 0.539586    
PRED_IND2   -0.882934    -0.383209   0.482871  -1.829 0.082427 .  
PRED_IND3   -1.451831    -0.553794   0.624294  -2.326 0.030675 *  
PRED_IND4   -1.151949    -0.439406   0.754374  -1.527 0.142415    
PRED_IND5   -1.098412    -0.459997   0.796382  -1.379 0.183047    
EXP_IND2     0.069167     0.027761   0.224152   0.309 0.760836    
EXP_IND3     0.027348     0.008247   0.293812   0.093 0.926767    
EXP_IND4    -0.385732    -0.116321   0.421313  -0.916 0.370812    
EXP_IND5     0.003063     0.001169   0.555897   0.006 0.995658    
DANG_IND2    0.703874     0.294771   0.457148   1.540 0.139305    
DANG_IND3    1.044139     0.398282   0.624193   1.673 0.109941    
DANG_IND4    1.275729     0.534255   0.811954   1.571 0.131827    
DANG_IND5    1.162274     0.386226   1.122186   1.036 0.312693    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3471 on 20 degrees of freedom
Multiple R-squared:  0.939, Adjusted R-squared:  0.884 
F-statistic: 17.09 on 18 and 20 DF,  p-value: 1.994e-08

Регрессия значима, но значимых коэффициентов на уровне \(\alpha = 0.05\) только два: для веса мозга и для фиктивной переменной индекса хищничества при значении \(3\).

14. Доверительные эллипсы

Посчитаем ковариационную и корреляционную матрицу для коэффициентов регрессии:

dummy.variable <- data.frame(NAME = df.log.out.na$NAME)

for (i in 2:5){
  dummy.variable[[paste("PRED_IND", i, sep = "_")]] = ifelse(df.log.out.na$PRED_IND == i, 1, 0)
}
for (i in 2:5){
  dummy.variable[[paste("EXP_IND", i, sep = "_")]] = ifelse(df.log.out.na$EXP_IND == i, 1, 0)
}
for (i in 2:5){
  dummy.variable[[paste("DANG_IND", i, sep = "_")]] = ifelse(df.log.out.na$DANG_IND == i, 1, 0)
}

n.df <- length(df.log.out.na$NAME)
sigma.df <- sum(model$residuals ** 2) / (n.df - model$rank)
covMatrix.df <- df.log.out.na |> dplyr::select(-NAME, -PRED_IND, -EXP_IND, -DANG_IND) |> cbind(dummy.variable |> dplyr::select(-NAME)) |> cov()
cov_b.df <- sigma.df / n.df * solve(covMatrix.df)
cov_beta.df <- sigma.df / n.df * solve(cov2cor(covMatrix.df))
cor_b.df <- cov2cor(cov_b.df)
corrplot(cor_b.df, method = "color")

Видим сильную корреляцию между между оценками коэффициентов регрессии между весом тела и мозга, между всеми видами сна, между индексами общей опасности и хищничества, а также между весом мозга и продолжительностью жизни.

Посмотрим на доверительный эллипсоид для значимых коэффициентов, обозначенных выше:

plot(ellipse::ellipse(cov_beta.df[c(2, 9), c(2, 9)], centre = model.beta$standardized.coefficients[c(3, 9)], level=0.95, npoints = 100), type = "l", asp = 1)
lines(x = c(-2, 2), y = c(2, -2), col = "red")

Он имеет плохой вид: либо оба имеет сильное влияние, либо оба — слабое. Причём индекс хищничества может быть совсем не значим, так как эллипс пересекает нулевое значение по \(y\).

15. Сильно коррелирующие признаки

Построим таблицу VIF:

ols_vif_tol(model)

Наблюдаем очень большую мультиколлинеарность почти по всем признакам.

Теперь таблица с частичными корреляциями:

ols_correlations(model)

                Correlations                 
--------------------------------------------
Variable     Zero Order    Partial     Part  
--------------------------------------------
BODY_WEI          0.842     -0.052    -0.013 
BRAIN_WE          0.933      0.703     0.244 
SLOWWAVE         -0.497      0.120     0.030 
PARADOX          -0.369      0.080     0.020 
SLEEP            -0.517     -0.062    -0.015 
GESTTIME          0.737     -0.138    -0.034 
PRED_IND2        -0.209     -0.378    -0.101 
PRED_IND3        -0.241     -0.461    -0.128 
PRED_IND4         0.075     -0.323    -0.084 
PRED_IND5         0.090     -0.295    -0.076 
EXP_IND2         -0.243      0.069     0.017 
EXP_IND3          0.220      0.021     0.005 
EXP_IND4          0.158     -0.201    -0.051 
EXP_IND5          0.455      0.001     0.000 
DANG_IND2        -0.215      0.326     0.085 
DANG_IND3        -0.399      0.350     0.092 
DANG_IND4         0.218      0.331     0.087 
DANG_IND5         0.307      0.226     0.057 
--------------------------------------------

Сильное влияние на зависимую переменную за исключением остальных оказывают, в перую очередь, вес мозга, а также индекс хищничества и общий индекс опасности (хотя он сильно зависит от индекса хищничества).

Попробуем убрать сильно коррелирующие. Уберём вес тела, общее время сна и общий индекс опасности:

model.minuscor <- lm(LIFESPAN ~ ., df.log.out.na |> dplyr::select(-NAME, -BODY_WEI, - SLEEP, -DANG_IND))
model.minuscor.beta <- lm.beta(model.minuscor)
summary(model.minuscor.beta)

Call:
lm(formula = LIFESPAN ~ ., data = dplyr::select(df.log.out.na, 
    -NAME, -BODY_WEI, -SLEEP, -DANG_IND))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.70225 -0.16574  0.03341  0.20078  0.47870 

Coefficients:
             Estimate Standardized Std. Error t value Pr(>|t|)    
(Intercept)  1.552951           NA   0.629838   2.466   0.0206 *  
BRAIN_WE     0.427375     1.015193   0.061132   6.991 2.01e-07 ***
SLOWWAVE     0.023400     0.080180   0.029182   0.802   0.4299    
PARADOX     -0.079520    -0.052792   0.153749  -0.517   0.6094    
GESTTIME    -0.087863    -0.092841   0.130827  -0.672   0.5078    
PRED_IND2   -0.371255    -0.161131   0.227439  -1.632   0.1147    
PRED_IND3   -0.534832    -0.204009   0.240393  -2.225   0.0350 *  
PRED_IND4   -0.062608    -0.023882   0.326654  -0.192   0.8495    
PRED_IND5   -0.053963    -0.022599   0.345540  -0.156   0.8771    
EXP_IND2     0.206412     0.082846   0.184211   1.121   0.2727    
EXP_IND3     0.121182     0.036543   0.258597   0.469   0.6433    
EXP_IND4    -0.189088    -0.057021   0.323711  -0.584   0.5642    
EXP_IND5    -0.020639    -0.007873   0.357456  -0.058   0.9544    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3648 on 26 degrees of freedom
Multiple R-squared:  0.9123,    Adjusted R-squared:  0.8719 
F-statistic: 22.55 on 12 and 26 DF,  p-value: 1.012e-10

По значимости имеем ту же ситуацию, однако значимость по p-value подросла, а \(R^2\) почти не изменился.

Мультиколлинеарность исчезла:

ols_vif_tol(model.minuscor)

16. Пошаговая регрессия

Backward

Построим пошаговую регрессию. Сначала рассмотрим убавление количества параметров:

st.bw.lm <- ols_step_backward_p(model.minuscor)
st.bw.lm

                            Stepwise Summary                             
-----------------------------------------------------------------------
Step    Variable       AIC       SBC       SBIC        R2       Adj. R2 
-----------------------------------------------------------------------
 0      Full Model    44.212    67.502    -71.965    0.91234    0.87188 
 1      PARADOX       42.612    64.238    -74.719    0.91144    0.87536 
 2      GESTTIME      40.951    60.914    -77.453    0.91066    0.87875 
 3      EXP_IND       35.968    49.276    -78.073    0.90348    0.88538 
 4      SLOWWAVE      34.481    46.126    -80.367    0.90220    0.88738 
-----------------------------------------------------------------------

Final Model Output 
------------------

                        Model Summary                          
--------------------------------------------------------------
R                       0.950       RMSE                0.315 
R-Squared               0.902       MSE                 0.117 
Adj. R-Squared          0.887       Coef. Var          14.153 
Pred R-Squared          0.864       AIC                34.481 
MAE                     0.260       SBC                46.126 
--------------------------------------------------------------
 RMSE: Root Mean Square Error 
 MSE: Mean Square Error 
 MAE: Mean Absolute Error 
 AIC: Akaike Information Criteria 
 SBC: Schwarz Bayesian Criteria 

                               ANOVA                                
-------------------------------------------------------------------
               Sum of                                              
              Squares        DF    Mean Square      F         Sig. 
-------------------------------------------------------------------
Regression     35.614         5          7.123    60.884    0.0000 
Residual        3.861        33          0.117                     
Total          39.475        38                                    
-------------------------------------------------------------------

                                  Parameter Estimates                                    
----------------------------------------------------------------------------------------
      model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
----------------------------------------------------------------------------------------
(Intercept)     1.534         0.175                  8.760    0.000     1.178     1.890 
   BRAIN_WE     0.375         0.024        0.890    15.606    0.000     0.326     0.424 
  PRED_IND2    -0.311         0.184       -0.135    -1.686    0.101    -0.685     0.064 
  PRED_IND3    -0.569         0.196       -0.217    -2.909    0.006    -0.968    -0.171 
  PRED_IND4    -0.150         0.193       -0.057    -0.777    0.443    -0.541     0.242 
  PRED_IND5    -0.114         0.183       -0.048    -0.625    0.536    -0.487     0.258 
----------------------------------------------------------------------------------------

Как это выглядит по шагам:

plot(st.bw.lm)

Какая же модель получается в итоге:

summary(st.bw.lm$model)

Call:
lm(formula = paste(response, "~", paste(c(include, cterms), collapse = " + ")), 
    data = l)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.7207 -0.2307 -0.0212  0.2462  0.5830 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.53397    0.17511   8.760 3.99e-10 ***
BRAIN_WE     0.37477    0.02401  15.606  < 2e-16 ***
PRED_IND2   -0.31054    0.18418  -1.686  0.10121    
PRED_IND3   -0.56931    0.19574  -2.909  0.00645 ** 
PRED_IND4   -0.14951    0.19252  -0.777  0.44291    
PRED_IND5   -0.11438    0.18295  -0.625  0.53612    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.342 on 33 degrees of freedom
Multiple R-squared:  0.9022,    Adjusted R-squared:  0.8874 
F-statistic: 60.88 on 5 and 33 DF,  p-value: 1.066e-15

Осталось \(2\) признака, удалились время беременности, быстрый сон, медленный сон и индекс воздействия сна.

Forward

Теперь будем добавлять признаки:

st.fw.lm <- ols_step_forward_p(model.minuscor)
st.fw.lm

                             Stepwise Summary                              
-------------------------------------------------------------------------
Step    Variable        AIC        SBC       SBIC        R2       Adj. R2 
-------------------------------------------------------------------------
 0      Base Model    115.149    118.476      1.181    0.00000    0.00000 
 1      BRAIN_WE       37.664     42.655    -73.233    0.86972    0.86620 
 2      PRED_IND       34.481     46.126    -80.767    0.90220    0.88738 
-------------------------------------------------------------------------

Final Model Output 
------------------

                        Model Summary                          
--------------------------------------------------------------
R                       0.950       RMSE                0.315 
R-Squared               0.902       MSE                 0.117 
Adj. R-Squared          0.887       Coef. Var          14.153 
Pred R-Squared          0.864       AIC                34.481 
MAE                     0.260       SBC                46.126 
--------------------------------------------------------------
 RMSE: Root Mean Square Error 
 MSE: Mean Square Error 
 MAE: Mean Absolute Error 
 AIC: Akaike Information Criteria 
 SBC: Schwarz Bayesian Criteria 

                               ANOVA                                
-------------------------------------------------------------------
               Sum of                                              
              Squares        DF    Mean Square      F         Sig. 
-------------------------------------------------------------------
Regression     35.614         5          7.123    60.884    0.0000 
Residual        3.861        33          0.117                     
Total          39.475        38                                    
-------------------------------------------------------------------

                                  Parameter Estimates                                    
----------------------------------------------------------------------------------------
      model      Beta    Std. Error    Std. Beta      t        Sig      lower     upper 
----------------------------------------------------------------------------------------
(Intercept)     1.534         0.175                  8.760    0.000     1.178     1.890 
   BRAIN_WE     0.375         0.024        0.890    15.606    0.000     0.326     0.424 
  PRED_IND2    -0.311         0.184       -0.135    -1.686    0.101    -0.685     0.064 
  PRED_IND3    -0.569         0.196       -0.217    -2.909    0.006    -0.968    -0.171 
  PRED_IND4    -0.150         0.193       -0.057    -0.777    0.443    -0.541     0.242 
  PRED_IND5    -0.114         0.183       -0.048    -0.625    0.536    -0.487     0.258 
----------------------------------------------------------------------------------------

Как это выглядело на графике:

plot(st.fw.lm)

Какую же модель получили:

summary(st.fw.lm$model)

Call:
lm(formula = paste(response, "~", paste(preds, collapse = " + ")), 
    data = l)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.7207 -0.2307 -0.0212  0.2462  0.5830 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.53397    0.17511   8.760 3.99e-10 ***
BRAIN_WE     0.37477    0.02401  15.606  < 2e-16 ***
PRED_IND2   -0.31054    0.18418  -1.686  0.10121    
PRED_IND3   -0.56931    0.19574  -2.909  0.00645 ** 
PRED_IND4   -0.14951    0.19252  -0.777  0.44291    
PRED_IND5   -0.11438    0.18295  -0.625  0.53612    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.342 on 33 degrees of freedom
Multiple R-squared:  0.9022,    Adjusted R-squared:  0.8874 
F-statistic: 60.88 on 5 and 33 DF,  p-value: 1.066e-15

Она аналогична той, которую получили от backward.

Однако только при одном значении индекса хищничества коэффициент значимый, поэтому вручную создадим фиктивные переменные и проведём опять автоматический отбор:

df.log.out.dummy <- df.log.out.na
for (i in 2:5){
  df.log.out.dummy[[paste("PRED_IND", i, sep = "_")]] <- ifelse(df.log.out.dummy$PRED_IND == i, 1, 0)
}
model.out.dummy <- lm(LIFESPAN ~ BRAIN_WE + PRED_IND_2 + PRED_IND_3 + PRED_IND_4 + PRED_IND_5, data = df.log.out.dummy)
st.bw.lm.dummy <- ols_step_backward_p(model.out.dummy)
plot(st.bw.lm.dummy)

summary(st.bw.lm.dummy$model)

Call:
lm(formula = paste(response, "~", paste(c(include, cterms), collapse = " + ")), 
    data = l)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.63197 -0.23515 -0.04439  0.22367  0.67817 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.42794    0.10776  13.251 3.35e-15 ***
BRAIN_WE     0.37812    0.02315  16.332  < 2e-16 ***
PRED_IND_2  -0.21197    0.13117  -1.616  0.11507    
PRED_IND_3  -0.47162    0.14733  -3.201  0.00291 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3354 on 35 degrees of freedom
Multiple R-squared:  0.9003,    Adjusted R-squared:  0.8917 
F-statistic: 105.3 on 3 and 35 DF,  p-value: < 2.2e-16

Осталось два значения индекса хищничества.

17. Остатки и Predicted vs Residuals

Нашли лучшую модель регрессии. Добавим те наблюдения, которые имели пропущенные значения только на убранных данных:

df.log.na.best <- df.log |> filter(!is.na(LIFESPAN))
for (i in 2:5){
  df.log.na.best[[paste("PRED_IND", i, sep = "_")]] <- ifelse(df.log.na.best$PRED_IND == i, 1, 0)
}
model.best <- lm(LIFESPAN ~ BRAIN_WE + PRED_IND_2 + PRED_IND_3, data = df.log.na.best)
model.beta.best <- lm.beta(model.best)
summary(model.beta.best)

Call:
lm(formula = LIFESPAN ~ BRAIN_WE + PRED_IND_2 + PRED_IND_3, data = df.log.na.best)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.12708 -0.35389 -0.02041  0.23742  1.77925 

Coefficients:
            Estimate Standardized Std. Error t value Pr(>|t|)    
(Intercept)  1.80029           NA    0.14546  12.377  < 2e-16 ***
BRAIN_WE     0.28962      0.76276    0.03015   9.605 2.79e-13 ***
PRED_IND_2  -0.22525     -0.10042    0.18422  -1.223  0.22676    
PRED_IND_3  -0.52686     -0.22082    0.19306  -2.729  0.00856 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5536 on 54 degrees of freedom
Multiple R-squared:  0.6738,    Adjusted R-squared:  0.6557 
F-statistic: 37.19 on 3 and 54 DF,  p-value: 3.589e-13

Ясно, что согласованность модели упала, так как мы построили её на большем наборе данных (плюс были добавлены выбросы, удаленные раньше).

Теперь посмотрим на нормальность остатков (хотим точность проверяемых критериев):

plot(model.best, which=2)

Видим, что справа нормальности не наблюдается, хотя на основном массиве данных нормальность присутствует.

Теперь посмотрим на график Predicted vs Residuals:

plot(model.best,which=1)

Видим достаточно равномерное распределение в прямоугольной области около нуля. То есть можно предположить адекватность выбранной линейной модели, а также гомогедостичность (то есть равенство дисперсий остатков).

График Residuls vs Deleted Residuals:

del_res <- sapply(1:58, function(n) df.log.na.best$LIFESPAN[n] - predict(lm(LIFESPAN ~ BRAIN_WE + PRED_IND_2 + PRED_IND_3, data = df.log.na.best[-n, ]), df.log.na.best[n, ]))
plot(x = model.best$residuals, y = del_res, xlab = "Rsiduals", ylab = "Deleted residuals")
lines(x = c(-3, 3), y = c(-3, 3))

abline(lm(del_res ~ model.best$residuals), col = "blue")
text(x = model.best$residuals, y = del_res, labels = rownames(df.log.na.best), cex = 0.6, pos = 4, col = "red")

Как и ожидалось, точки расположились под чуть большим наклоном, нежели \(y = x\). Больших отклонений (по линии регресии на этих точках) не видно, то есть сказать что-то о выбросах нельзя.

18. Выбросы по Махаланобису

df.log.na.best.x <- df.log.na.best |> dplyr::select(BRAIN_WE, PRED_IND_2, PRED_IND_3)

cov_matrix <- cov(df.log.na.best.x)

mahalanobis_distance <- mahalanobis(df.log.na.best.x, center = colMeans(df.log.na.best.x), cov = cov_matrix)

plot(mahalanobis_distance, 
     main = "Mahalanobis Distance Plot",
     xlab = "Observation",
     ylab = "Mahalanobis Distance",
     type = "h", 
     col = "darkblue") 

plot(sort(mahalanobis_distance), 
     main = "Mahalanobis Distance Plot",
     xlab = "Observation",
     ylab = "Mahalanobis Distance",
     type = "h", 
     col = "darkblue") 

По Махаланобису виден скачок для двух наблюдений: \(1\) и \(4\). Это два слона:

NAME	BODY_WEI	BRAIN_WE	SLOWWAVE	PARADOX	SLEEP	LIFESPAN	GESTTIME	PRED_IND	EXP_IND	DANG_IND	PRED_IND_2	PRED_IND_3	PRED_IND_4	PRED_IND_5
African elephant	8.802974	8.650324	NA	NA	3.3	3.653252	6.46925	3	5	3	0	1	0	0
Asian elephant	7.842671	8.434463	2.1	0.5877867	3.9	4.234107	6.43615	3	5	4	0	1	0	0

19. Выбросы по Куку

cooks_distance <- cooks.distance(model.best)

plot(cooks_distance, 
     col = "darkblue",
     type = "h", 
     main = "Cook's Distance",
     xlab = "Observation",
     ylab = "Cook's Distance")

plot(sort(cooks_distance), 
     col = "darkblue",
     type = "h", 
     main = "Cook's Distance",
     xlab = "Observation",
     ylab = "Cook's Distance")

По скачку значения расстояния Кука можно выделить три наблюдения: 6, 14 и 31. Это те самые летучие мыши и ещё ехидна:

df.log.na.best[c(6, 14, 31), ] |> print_df()

NAME	BODY_WEI	BRAIN_WE	SLOWWAVE	PARADOX	SLEEP	LIFESPAN	GESTTIME	PRED_IND	EXP_IND	DANG_IND	PRED_IND_2	PRED_IND_3	PRED_IND_4	PRED_IND_5
Big brown bat	-3.772261	-1.203973	15.8	1.3609766	19.7	2.944439	3.555348	1	1	1	0	0	0	0
Echidna	1.098612	3.218876	8.6	NA	8.6	3.912023	3.332205	2	2	2	1	0	0	0
Little brown bat	-4.605170	-1.386294	17.9	0.6931472	19.9	3.178054	3.912023	1	1	1	0	0	0	0

20. Studentized Residuals vs Leverage Plot

Наконец посмотрим на соотношение рычага и остатков:

ols_plot_resid_lev(model.best)

Здесь явных выбросов нет: по остаткам, за пределами доверительного интервала небольшое количество наблюдений, по рычагу два кандидата на выбросы — слоны, а наблюдений, которые бы были и по рычагу и по остаткам “выбросами”, нет.

21. Исключаем выбросы

Найденные выбросы имеют определённый смысл — большие и средне атакуемые слоны, маленькие, но долгоживущие летучие мыши и такая же ехидна. Слонов исключать не будем, так как они просто большие, поэтому Махаланобис их выделяет, однако они хорошо согласуются с моделью (видно по картинке, например). Исключим выбросы:

model.best.out <- lm(LIFESPAN ~ BRAIN_WE + PRED_IND_2 + PRED_IND_3, data = df.log.na.best[c(-6, -14, -31),])
model.best.out.beta <- lm.beta(model.best.out)
summary(model.best.out.beta)

Call:
lm(formula = LIFESPAN ~ BRAIN_WE + PRED_IND_2 + PRED_IND_3, data = df.log.na.best[c(-6, 
    -14, -31), ])

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0491 -0.2596  0.0245  0.2405  1.0616 

Coefficients:
            Estimate Standardized Std. Error t value Pr(>|t|)    
(Intercept)  1.47436           NA    0.11487  12.835   <2e-16 ***
BRAIN_WE     0.34604      0.87394    0.02316  14.941   <2e-16 ***
PRED_IND_2  -0.15023     -0.06611    0.13724  -1.095    0.279    
PRED_IND_3  -0.36992     -0.15766    0.13833  -2.674    0.010 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3905 on 51 degrees of freedom
Multiple R-squared:  0.8394,    Adjusted R-squared:   0.83 
F-statistic: 88.87 on 3 and 51 DF,  p-value: < 2.2e-16

Судя по \(R^2\) и его исправленному варианту, а также по p-value значимости регрессии, получили более хорошую модель.

22. Предсказания

Предскажем время жизниобыкновенного ежа. Вес его мозга — \(3.3\) г, имеет защиту от хищников, но крупные хищные птицы успешно их атакуют, поэтому индекс хищничества — примерно \(3\). В неволе ежи живут \(4\)-\(6\) лет. Пусть будет \(5\):

data.hedgehog <- data.frame(NAME = "European hedgehog", BRAIN_WE = log(3.3), LIFESPAN = log(5), PRED_IND_2 = 0, PRED_IND_3 = 1)
data.hedgehog |> print_df()

NAME	BRAIN_WE	LIFESPAN	PRED_IND_2	PRED_IND_3
European hedgehog	1.193922	1.609438	0	1

Предскажем его продолжительность жизни:

pred.hedgehog.pred <- predict(model.best.out, newdata = data.hedgehog, interval = "predict")
pred.hedgehog.conf <- predict(model.best.out, newdata = data.hedgehog, interval = "confidence")

exp(pred.hedgehog.pred[, 1])

[1] 4.561225

Получили цифры близкие к действительности. Посмотрим на предсказательный интервал:

print("Нижняя граница:")

[1] "Нижняя граница:"

exp(pred.hedgehog.pred[, 2])

[1] 2.002536

print("Верхняя граница:")

[1] "Верхняя граница:"

exp(pred.hedgehog.pred[, 3])

[1] 10.38922

Ежи в неволе могут жить до 10 лет.

А теперь на доверительный:

print("Нижняя граница:")

[1] "Нижняя граница:"

exp(pred.hedgehog.conf[, 2])

[1] 3.549484

print("Верхняя граница:")

[1] "Верхняя граница:"

exp(pred.hedgehog.conf[, 3])

[1] 5.861353